//******************************************************************************
// mumupi0pi0_h5.cc - A BASF module for reading MDST files and extracting 
// quantities of interest for events consistent with a mu+mu-pi0pi0 final state.
//
// Date: June 2012
// Author: craigb (craig dot bookwalter at pnnl dot gov)
//******************************************************************************

// Standard headers
#include <algorithm>
#include <exception>
#include <iostream>
#include <fstream>
#include <set>
#include <stdexcept>
#include <string>
#include <sstream>
#include <utility>
#include <vector>

// Belle-specific headers
#include "panther/panther.h"
#include MDST_H
#include BELLETDF_H
#include HEPEVT_H

#include "particle/Particle.h"
#include "event/BelleEvent.h"
#include "tuple/BelleTupleManager.h"
#include "basf/module.h"
#include "basf/module_descr.h"
#include "benergy/BeamEnergy.h"
#include "eid/eid.h"


// Headers written by me
#include <pytpp.hpp>

#include "rows/ypipi_pi0pi0row.hpp"
#include "ypipi_utils.hpp"

using std::string;
using std::vector;
using std::cout;
using std::endl;

static const double RADS_TO_DEG = 360./(2*3.14159265359);

#ifdef BELLE_NAMESPACE
namespace Belle {
#endif //BELLE_NAMESPACE

  // All BASF modules inherit from the "Module" class, defined in `basf/module.h'
  // They are written as a single class. The BASF executable loads this class
  // (compiled into a shared object) dynamically. The module constructor, `init' 
  // and `hist_def' methods are all called before reading events from the event 
  // source. The first event is read and BASF calls `begin_run()'. Every time 
  // the experimental run number (contained in the datastream) changes, `end_run()' 
  // will be called and then `begin_run()' again.
  //
  // The `event()' function does the event-by-event processing; in the case of 
  // this module, saving some quantities of interest to an HDF5 dataset for 
  // later study. 
  //
  // When the event source is exhausted, the term() function is called, followed 
  // by the class destructor.
  //----------------------------------------------------------------------------
  class pi0pi0inclusive_h5: public Module {
  public:
    pi0pi0inclusive_h5();
    virtual ~pi0pi0inclusive_h5();
    void init(int *);
    void term();
    void disp_stat(const char*);
    void hist_def();
    void event(BelleEvent*, int*);
    void hbook_event(BelleEvent*, int*);
    void h5event(BelleEvent*, int*);
    void begin_run(BelleEvent*, int*);
    void end_run(BelleEvent*, int*);

    char outpath_stub[4096];
    char inpath[4096]; 
    char nttype[256];
    char output_format[256];

    BelleTuple *_bt;

  private:
    float calc_R2();

    bool _write_h5;
    bool _write_hbook;

    int _nevents_tot;
    int _event_yield;

    pyt::table<ypipi_pi0pi0row_t>* _nt; 
    ypipi_pi0pi0row_t* _row;
  };


  //----------------------------------------------------------------------------
  // This function registers the module with the BASF executeable.
  // Any parameters defined in the control script (`pi0pi0inclusive_h5.basf') are
  // declared and matched with public member variables of the module class.
  //----------------------------------------------------------------------------
  extern "C" Module_descr *mdcl_pi0pi0inclusive_h5()
  {

    pi0pi0inclusive_h5 *tmp = new pi0pi0inclusive_h5();
    Module_descr *dscr = new Module_descr ( "pi0pi0inclusive_h5", tmp);
    dscr->define_param("outpath_stub", "Path stub to output file", "S", 4096, tmp->outpath_stub);
    dscr->define_param("inpath", "Input file", "S", 4096, tmp->inpath);
    dscr->define_param("nttype", "N-tuple type", "S", 256, tmp->nttype);
    dscr->define_param("output_format", "Output format", "S", 256, tmp->output_format);

    BeamEnergy::define_global(dscr);

    return dscr;
  }
  
  //----------------------------------------------------------------------------
  // Default constructor
  //----------------------------------------------------------------------------
  pi0pi0inclusive_h5::pi0pi0inclusive_h5() : 
  _nevents_tot(0),
  _event_yield(0), 
    _row(new ypipi_pi0pi0row_t),
    _bt(NULL),
    _nt(NULL),
    _write_h5(false),
    _write_hbook(false)
  {
    // Initialize member vars
  }

  //----------------------------------------------------------------------------
  // Default destructor
  //----------------------------------------------------------------------------  
  pi0pi0inclusive_h5::~pi0pi0inclusive_h5() {
    if (_write_h5) {
      delete _nt;
      delete _row;
    }
  }

  //----------------------------------------------------------------------------
  // init() I leave empty
  //----------------------------------------------------------------------------    
  void pi0pi0inclusive_h5::init(int* status) {
    if (strstr(output_format, "hbook"))
      _write_hbook = true;

    if (strstr(output_format, "h5"))
      _write_h5 = true;

    string path_stub(outpath_stub);

    if (_write_h5) {
      string outpath = path_stub + ".h5";
      
      _nt = pyt::open_table<ypipi_pi0pi0row_t>(outpath, "nt", "w");
    
      _nt->table_setattr("mdst_path", inpath);
      _nt->table_setattr("nttype", nttype);
      _nt->table_setattr("module", "ypipi");
    }
  }
  
  //----------------------------------------------------------------------------
  // in hist_def(), I define one or more output files based on the values of the 
  // parameters passed in from the control scripts.
  //----------------------------------------------------------------------------  
  void pi0pi0inclusive_h5::hist_def() {
    extern BelleTupleManager* BASF_Histogram;

    string path_stub(outpath_stub);

    if (_write_hbook) {
      _bt = BASF_Histogram->ntuple("_bt", "expno runno eventno ebeam");
    }

    if ((_nt == NULL) && (_bt == NULL)) {
      std::stringstream ss; 
      ss << "Output format `" << output_format << "' does not contain `hbook' or `h5' -- no output will be written." 
	 << endl << "Aborting...";

      throw std::invalid_argument(ss.str().c_str());
    }

    cout << "CRAIGB : Analyzing file " << inpath << endl;

  }


  //----------------------------------------------------------------------------
  // begin_run() just contains some bookkeeping to make sure the numbers I save
  // are well-calibrated.
  //----------------------------------------------------------------------------    
  void pi0pi0inclusive_h5::begin_run(BelleEvent* be, int* status) {
    eid::init_data();
    BeamEnergy::begin_run();
  }
 

  void pi0pi0inclusive_h5::hbook_event(BelleEvent* evptr, int* status) {
    Belle_event_Manager& m_event = Belle_event_Manager::get_manager();
    
    _bt->column("expno", m_event[0].ExpNo());
    _bt->column("runno", m_event[0].RunNo());
    _bt->column("eventno", m_event[0].EvtNo());
    _bt->column("ebeam", BeamEnergy::E_beam_corr());

    _bt->dumpData();
    _bt->clearData();
  }
    
  //----------------------------------------------------------------------------
  // hbook_event() is called if the `hbook' output format is requested from the
  // control script; it processes a single event and saves quantities of 
  // interest to the output file.
  //----------------------------------------------------------------------------  
  // void pi0pi0inclusive_h5::hbook_event(BelleEvent* evptr, int* status) {
  //   Belle_event_Manager& m_event = Belle_event_Manager::get_manager();
 
  //   Hep3Vector cm_boost = BeamEnergy::CMBoost();
  //   HepLorentzVector initial(0, 0, 0, 10.876);

  //   Mdst_pi0_Manager& m_pi0 = Mdst_pi0_Manager::get_manager();
  //   Mdst_ecl_aux_Manager& m_ecl_aux = Mdst_ecl_aux_Manager::get_manager();

  //   vector<Mdst_pi0> pi0_candidates;

  //   for (Mdst_pi0_Manager::iterator it = m_pi0.begin();
  // 	 it != m_pi0.end(); it++) 
  //     {
  // 	// if (fabs(0.135 - it->mass()) >= 0.0054*3) 
  // 	//   continue;

  // 	pi0_candidates.push_back(*it);
  //     }

  //   if (pi0_candidates.size() < 2)
  //     return;
    
  //   vector<std::pair<Mdst_pi0, Mdst_pi0> > pi0pi0_combinations;

  //   for (unsigned int i=0; i < pi0_candidates.size(); i++) 
  //     for (unsigned int j=i+1; j < pi0_candidates.size(); j++) 
  // 	pi0pi0_combinations.push_back(std::pair<Mdst_pi0, Mdst_pi0>(pi0_candidates[i], pi0_candidates[j]));

  //   std::sort(pi0pi0_combinations.begin(), pi0pi0_combinations.end(), pi0pi0_chisq_cmp);
    
  //   for (vector<std::pair<Mdst_pi0, Mdst_pi0> >::iterator it = pi0pi0_combinations.begin();
  // 	 it != pi0pi0_combinations.end(); it++)
  //     {
  // 	Mdst_pi0& pi01_mdst = it->first;
  // 	Mdst_pi0& pi02_mdst = it->second;
	  
  // 	if (it->first.energy() < it->second.energy()) {
  // 	  Mdst_pi0& pi01_mdst = it->second;
  // 	  Mdst_pi0& pi02_mdst = it->first;
  // 	}

  // 	std::set<Panther_ID> gamma_IDs;

  // 	gamma_IDs.insert(pi01_mdst.gamma(0).get_ID());
  // 	gamma_IDs.insert(pi01_mdst.gamma(1).get_ID());
  // 	gamma_IDs.insert(pi02_mdst.gamma(0).get_ID());
  // 	gamma_IDs.insert(pi02_mdst.gamma(1).get_ID());
	
  // 	if (gamma_IDs.size() != 4)
  // 	  continue;

  // 	_bt->column("expno", m_event[0].ExpNo());
  // 	_bt->column("runno", m_event[0].RunNo());
  // 	_bt->column("eventno", m_event[0].EvtNo());
  // 	_bt->column("ebeam", BeamEnergy::E_beam_corr());
	
  // 	_bt->column("r2", calc_R2());
	
  // 	_bt->column("pi01chi2", pi01_mdst.chisq());
  // 	_bt->column("pi02chi2", pi02_mdst.chisq());
  // 	_bt->column("pi01m", pi01_mdst.mass());
  // 	_bt->column("pi02m", pi02_mdst.mass());

  // 	HepLorentzVector pi01(pi01_mdst.px(), pi01_mdst.py(), pi01_mdst.pz(), pi01_mdst.energy());
  // 	HepLorentzVector pi02(pi02_mdst.px(), pi02_mdst.py(), pi02_mdst.pz(), pi02_mdst.energy());

  // 	pi01.boost(-cm_boost);
  // 	pi02.boost(-cm_boost);

  // 	_bt->column("pi01cmpx", pi01.px());
  // 	_bt->column("pi01cmpy", pi01.py());
  // 	_bt->column("pi01cmpz", pi01.pz());
  // 	_bt->column("pi01cme", pi01.e());

  // 	_bt->column("pi02cmpx", pi02.px());
  // 	_bt->column("pi02cmpy", pi02.py());
  // 	_bt->column("pi02cmpz", pi02.pz());
  // 	_bt->column("pi02cme", pi02.e());

  // 	_bt->column("mmpi0pi0", (initial - pi01 - pi02).m());
  // 	_bt->column("mm2pi01", (initial - pi01).m2());
  // 	_bt->column("mm2pi02", (initial - pi02).m2());
  // 	_bt->column("mpi0pi0", (pi01 + pi02).m());
	
  // 	Mdst_ecl& pi01_g1ecl = pi01_mdst.gamma(0).ecl();
  // 	Mdst_ecl& pi01_g2ecl = pi01_mdst.gamma(1).ecl();

  // 	Hep3Vector pi01_g1(pi01_mdst.gamma(0).px(), 
  // 			   pi01_mdst.gamma(0).py(), 
  // 			   pi01_mdst.gamma(0).pz());

  // 	Hep3Vector pi01_g2(pi01_mdst.gamma(1).px(), 
  // 			   pi01_mdst.gamma(1).py(), 
  // 			   pi01_mdst.gamma(1).pz());

  // 	_bt->column("pi01g1th", pi01_g1.theta() * RADS_TO_DEG);
  // 	_bt->column("pi01g2th", pi01_g2.theta() * RADS_TO_DEG);
	  
  // 	if (pi01_g1ecl.energy() < pi01_g2ecl.energy()) {
  // 	  Mdst_ecl tmp = pi01_g2ecl;
  // 	  pi01_g2ecl = pi01_g1ecl;
  // 	  pi01_g1ecl = tmp;
  // 	}	 

  // 	_bt->column("pi01g1e", pi01_g1ecl.energy());
  // 	_bt->column("pi01g2e", pi01_g2ecl.energy());
	
  // 	Mdst_ecl_aux& pi01_eclaux_g1 = m_ecl_aux(pi01_g1ecl.get_ID());
  // 	Mdst_ecl_aux& pi01_eclaux_g2 = m_ecl_aux(pi01_g2ecl.get_ID());

  // 	_bt->column("pi01g1e9", pi01_eclaux_g1.e9oe25());
  // 	_bt->column("pi01g2e9", pi01_eclaux_g2.e9oe25());

  // 	Hep3Vector pi02_g1(pi02_mdst.gamma(0).px(), 
  // 			   pi02_mdst.gamma(0).py(), 
  // 			   pi02_mdst.gamma(0).pz());

  // 	Hep3Vector pi02_g2(pi02_mdst.gamma(1).px(), 
  // 			   pi02_mdst.gamma(1).py(), 
  // 			   pi02_mdst.gamma(1).pz());

  // 	_bt->column("pi02g1th", pi02_g1.theta() * RADS_TO_DEG);
  // 	_bt->column("pi02g2th", pi02_g2.theta() * RADS_TO_DEG);

  // 	Mdst_ecl& pi02_g1ecl = pi02_mdst.gamma(0).ecl();
  // 	Mdst_ecl& pi02_g2ecl = pi02_mdst.gamma(1).ecl();

  // 	if (pi02_g1ecl.energy() < pi02_g2ecl.energy()) {
  // 	  Mdst_ecl tmp = pi02_g2ecl;
  // 	  pi02_g2ecl = pi02_g1ecl;
  // 	  pi02_g1ecl = tmp;
  // 	}	 

  // 	_bt->column("pi02g1e", pi02_g1ecl.energy());
  // 	_bt->column("pi02g2e", pi02_g2ecl.energy());
	
  // 	Mdst_ecl_aux& pi02_eclaux_g1 = m_ecl_aux(pi02_g1ecl.get_ID());
  // 	Mdst_ecl_aux& pi02_eclaux_g2 = m_ecl_aux(pi02_g2ecl.get_ID());
  // 	_bt->column("pi02g1e9", pi02_eclaux_g1.e9oe25());
  // 	_bt->column("pi02g2e9", pi02_eclaux_g2.e9oe25());

  // 	_bt->column("missinge", pi01.e() + pi02.e() + mup.e() + mum.e() - 2*BeamEnergy::E_beam_corr());

  // 	_bt->dumpData();

  // 	if (strcmp(this->nttype, "best_chi2") == 0)
  // 	  break;
  //     }

  //   _bt->clearData();
  //   _event_yield++;

  // }

  //----------------------------------------------------------------------------
  // h5event() is called if the `h5' output format is requested from the
  // control script; again, it processes a single event and saves quantities of 
  // interest to the output file.
  //----------------------------------------------------------------------------  
  void pi0pi0inclusive_h5::h5event(BelleEvent* evptr, int* status) {
    Belle_event_Manager& m_event = Belle_event_Manager::get_manager();

    _row->expno   = m_event[0].ExpNo();
    _row->runno   = m_event[0].RunNo();
    _row->eventno = m_event[0].EvtNo();
    _row->beam_e  = BeamEnergy::E_beam_corr();

    Hep3Vector boost_cm = BeamEnergy::CMBoost();
    HepLorentzVector initial(0, 0, 0, 10.876);

    vector<Particle> v_pi0;
  
    _row->r2 = calc_R2();
    
    Mdst_pi0_Manager& m_pi0 = Mdst_pi0_Manager::get_manager();
    Mdst_ecl_aux_Manager& m_ecl_aux = Mdst_ecl_aux_Manager::get_manager();

    vector<Mdst_pi0> pi0_candidates;

    for (Mdst_pi0_Manager::iterator it = m_pi0.begin();
    	 it != m_pi0.end(); it++) 
      {
    	// if (fabs(0.135 - it->mass()) >= 0.0054*3) 
    	//   continue;

    	pi0_candidates.push_back(*it);
      }

    if (pi0_candidates.size() < 2) 
      return;
    
    vector<std::pair<Mdst_pi0, Mdst_pi0> > pi0pi0_combinations;

    for (unsigned int i=0; i < pi0_candidates.size(); i++) 
      for (unsigned int j=i+1; j < pi0_candidates.size(); j++) 
    	pi0pi0_combinations.push_back(std::pair<Mdst_pi0, Mdst_pi0>(pi0_candidates[i], pi0_candidates[j]));

    std::sort(pi0pi0_combinations.begin(), pi0pi0_combinations.end(), pi0pi0_chisq_cmp);
    
    for (vector<std::pair<Mdst_pi0, Mdst_pi0> >::iterator it = pi0pi0_combinations.begin();
    	 it != pi0pi0_combinations.end(); it++)
      {

    	Mdst_pi0& pi01_mdst = it->first;
    	Mdst_pi0& pi02_mdst = it->second;
	  
    	if (it->first.energy() < it->second.energy()) {
    	  Mdst_pi0& pi01_mdst = it->second;
    	  Mdst_pi0& pi02_mdst = it->first;
    	}

    	std::set<Panther_ID> gamma_IDs;

    	gamma_IDs.insert(pi01_mdst.gamma(0).get_ID());
    	gamma_IDs.insert(pi01_mdst.gamma(1).get_ID());
    	gamma_IDs.insert(pi02_mdst.gamma(0).get_ID());
    	gamma_IDs.insert(pi02_mdst.gamma(1).get_ID());
	
    	if (gamma_IDs.size() != 4)
    	  continue;

    	_row->pi01_chisq = pi01_mdst.chisq();
    	_row->pi02_chisq = pi02_mdst.chisq();
    	_row->pi01_m     = pi01_mdst.mass();
    	_row->pi02_m     = pi02_mdst.mass();

    	HepLorentzVector pi01(pi01_mdst.px(), pi01_mdst.py(), pi01_mdst.pz(), pi01_mdst.energy());
    	HepLorentzVector pi02(pi02_mdst.px(), pi02_mdst.py(), pi02_mdst.pz(), pi02_mdst.energy());

    	pi01.boost(-boost_cm);
    	pi02.boost(-boost_cm);

    	_row->pi01_px_cm = pi01.px();
    	_row->pi01_py_cm = pi01.py();
    	_row->pi01_pz_cm = pi01.pz();
    	_row->pi01_e_cm  = pi01.e();

    	_row->pi02_px_cm = pi02.px();
    	_row->pi02_py_cm = pi02.py();
    	_row->pi02_pz_cm = pi02.pz();
    	_row->pi02_e_cm  = pi02.e();

	HepLorentzVector msng = initial - pi01 - pi02;
	_row->missing_px_cm = msng.px();
	_row->missing_py_cm = msng.py();
	_row->missing_pz_cm = msng.pz();
	_row->missing_e_cm = msng.e();

    	_row->mm_off_pi0pi0 = msng.m();
    	_row->m2_off_pi01 = (initial - pi01).m2();
    	_row->m2_off_pi02 = (initial - pi02).m2();
    	_row->m_pi0pi0 = (pi01 + pi02).m();
	
    	Mdst_ecl& pi01_g1ecl = pi01_mdst.gamma(0).ecl();
    	Mdst_ecl& pi01_g2ecl = pi01_mdst.gamma(1).ecl();

    	Hep3Vector pi01_g1(pi01_mdst.gamma(0).px(), 
    			   pi01_mdst.gamma(0).py(), 
    			   pi01_mdst.gamma(0).pz());

    	Hep3Vector pi01_g2(pi01_mdst.gamma(1).px(), 
    			   pi01_mdst.gamma(1).py(), 
    			   pi01_mdst.gamma(1).pz());

    	_row->pi01_g1theta_lab = pi01_g1.theta() * RADS_TO_DEG;
    	_row->pi01_g2theta_lab = pi01_g2.theta() * RADS_TO_DEG;
	  
    	if (pi01_g1ecl.energy() < pi01_g2ecl.energy()) {
    	  Mdst_ecl tmp = pi01_g2ecl;
    	  pi01_g2ecl = pi01_g1ecl;
    	  pi01_g1ecl = tmp;
    	}	 

    	_row->pi01_g1e = pi01_g1ecl.energy();
    	_row->pi01_g2e = pi01_g2ecl.energy();
	
    	Mdst_ecl_aux& pi01_eclaux_g1 = m_ecl_aux(pi01_g1ecl.get_ID());
    	Mdst_ecl_aux& pi01_eclaux_g2 = m_ecl_aux(pi01_g2ecl.get_ID());

    	_row->pi01_g1e9oe25 = pi01_eclaux_g1.e9oe25();
    	_row->pi01_g2e9oe25 = pi01_eclaux_g2.e9oe25();

    	Hep3Vector pi02_g1(pi02_mdst.gamma(0).px(), 
    			   pi02_mdst.gamma(0).py(), 
    			   pi02_mdst.gamma(0).pz());

    	Hep3Vector pi02_g2(pi02_mdst.gamma(1).px(), 
    			   pi02_mdst.gamma(1).py(), 
    			   pi02_mdst.gamma(1).pz());

    	_row->pi02_g1theta_lab = pi02_g1.theta() * RADS_TO_DEG;
    	_row->pi02_g2theta_lab = pi02_g2.theta() * RADS_TO_DEG;

    	Mdst_ecl& pi02_g1ecl = pi02_mdst.gamma(0).ecl();
    	Mdst_ecl& pi02_g2ecl = pi02_mdst.gamma(1).ecl();

    	if (pi02_g1ecl.energy() < pi02_g2ecl.energy()) {
    	  Mdst_ecl tmp = pi02_g2ecl;
    	  pi02_g2ecl = pi02_g1ecl;
    	  pi02_g1ecl = tmp;
    	}	 

    	_row->pi02_g1e = pi02_g1ecl.energy();
    	_row->pi02_g2e = pi02_g2ecl.energy();
	
    	Mdst_ecl_aux& pi02_eclaux_g1 = m_ecl_aux(pi02_g1ecl.get_ID());
    	Mdst_ecl_aux& pi02_eclaux_g2 = m_ecl_aux(pi02_g2ecl.get_ID());
    	_row->pi02_g1e9oe25 = pi02_eclaux_g1.e9oe25();
    	_row->pi02_g2e9oe25 = pi02_eclaux_g2.e9oe25();
		
    	_nt->append(*_row);	

    	if (strcmp(this->nttype, "best_chi2") == 0)
    	  break;
      }
     
    _event_yield++;       
  }

  //----------------------------------------------------------------------------
  // event() is called for each event read from the event source; in this module
  // it simply dispatches to different functions for saving to different output
  // formats.
  //----------------------------------------------------------------------------  
  void pi0pi0inclusive_h5::event(BelleEvent *evptr, int *status) {
    _nevents_tot++;

    if (_write_hbook)
      hbook_event(evptr, status);
    
    if (_write_h5) 
      h5event(evptr, status);
    
  }


  //----------------------------------------------------------------------------
  // calc_R2() is a function that calculates a physics quantity (normalized Fox-
  // Wolfram second moment, if you're curious).
  //----------------------------------------------------------------------------    
  float pi0pi0inclusive_h5::calc_R2() {
    Mdst_charged_Manager& m_charged = Mdst_charged_Manager::get_manager();
    Mdst_gamma_Manager& m_gamma = Mdst_gamma_Manager::get_manager();
    
    
    int n_trk = m_charged.count() + m_gamma.count();
    float H0, R1, R2, R3, R4;
    float p_trk[n_trk*3];

    float beam_e = BeamEnergy::E_beam_corr();

    int trkno = 0;
    for (Mdst_charged_Manager::iterator it = m_charged.begin(); 
	 it != m_charged.end(); it++, trkno++)
      {
	p_trk[3*trkno]     = it->px();
	p_trk[3*trkno + 1] = it->py();
	p_trk[3*trkno + 2] = it->pz();
      }

    for (Mdst_gamma_Manager::iterator it = m_gamma.begin();
	 it != m_gamma.end(); it++, trkno++)
      {
	p_trk[3*trkno]     = it->px();
	p_trk[3*trkno + 1] = it->py();
	p_trk[3*trkno + 2] = it->pz();
      }	

    fwjet2(n_trk, p_trk, beam_e, &H0, &R1, &R2, &R3, &R4);

    return R2; 
  }

  
  //----------------------------------------------------------------------------
  // end_run() handles any bookkeeping necessary at the end of an experimental
  // run--in this case it's empty.
  //----------------------------------------------------------------------------    
  void pi0pi0inclusive_h5::end_run(BelleEvent* be, int* status) { 
    // Yessir
  }

  //----------------------------------------------------------------------------
  // term() is used to write out any output waiting to be flushed
  //----------------------------------------------------------------------------    
  void pi0pi0inclusive_h5::term() {
    // if ((_write_hbook) && (_write_h5))
    //   _event_yield /= 2;

    cout << "CRAIGB : " << _nevents_tot << " events processed, " <<  _event_yield << " pi0pi0 events found" << endl;
    // End of analysis, close ntuples
  }


  //----------------------------------------------------------------------------
  // disp_stat() would display some statistics about the events I saved or threw
  // away if I was motivated enough to keep track of such things.
  //----------------------------------------------------------------------------    
  void pi0pi0inclusive_h5::disp_stat(const char* str) {
    // Who knows
  }

#ifdef BELLE_NAMESPACE
}
#endif
